library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(knitr)
library(ggplot2)
library(dplyr)
library(lubridate)
library(MMWRweek)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
library(readxl)
library(excessmort)
library(plotly)
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
library(patchwork)
library(scales)
## 
## Attaching package: 'scales'
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## The following object is masked from 'package:readr':
## 
##     col_factor
library(cowplot)
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:patchwork':
## 
##     align_plots
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
source("~/Desktop/Thesis Work/Hospitalization/11-26 data/Deficit_Model_1.0.R")
allcause_state_byweek <- read_excel("allcause_state_byweek_suppressed_25mar21.xlsx")
icd10categories <- read_excel("icd10categories_updated_321.xlsx")

#weekly state data
all_cause_state_byweek <- allcause_state_byweek %>%
  rename(weekly_count_statelevel = count)

all_cause_state_by_week <- all_cause_state_byweek %>%
  separate(mmwr_week, into = c("MMWRyear", "MMWRweek"), sep = "-", convert = TRUE) %>%
  mutate(
    date = MMWRweek2Date(MMWRyear = MMWRyear, MMWRweek = MMWRweek, MMWRday = 1),
    statecount_numeric = as.numeric(ifelse(weekly_count_statelevel == "<5", NA, weekly_count_statelevel)),
    statecount_suppressed = weekly_count_statelevel == "<5",
    category = as.character(category)
  )

all_cause_state_by_week %>%
  group_by(category) %>%
  reframe(missing_total = sum(statecount_suppressed), 
            count_total = sum(statecount_numeric, na.rm = TRUE))
## # A tibble: 10 × 3
##    category               missing_total count_total
##    <chr>                          <int>       <dbl>
##  1 cardiovascular                     0      373147
##  2 covid                              1       58262
##  3 digestive                          0      182701
##  4 infectious                         0      112279
##  5 infectious respiratory             0       29816
##  6 injury                             0      331205
##  7 other                              0     1740863
##  8 pneumonia                          0       97954
##  9 renal                              0      102199
## 10 respiratory                        0      164424
all_cause_state_by_week
## # A tibble: 3,140 × 7
##    MMWRyear MMWRweek category               weekly_count_statelevel date      
##       <int>    <int> <chr>                  <chr>                   <date>    
##  1     2019        1 cardiovascular         843                     2018-12-30
##  2     2019        1 digestive              441                     2018-12-30
##  3     2019        1 infectious             299                     2018-12-30
##  4     2019        1 infectious respiratory 240                     2018-12-30
##  5     2019        1 injury                 691                     2018-12-30
##  6     2019        1 other                  4090                    2018-12-30
##  7     2019        1 pneumonia              427                     2018-12-30
##  8     2019        1 renal                  220                     2018-12-30
##  9     2019        1 respiratory            522                     2018-12-30
## 10     2019        1 covid                  0                       2018-12-30
## # ℹ 3,130 more rows
## # ℹ 2 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>
str(icd10categories)
## tibble [869 × 5] (S3: tbl_df/tbl/data.frame)
##  $ icd10code          : chr [1:869] "A00" "A01" "A02" "A03" ...
##  $ category_original  : chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
##  $ category_1113update: chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
##  $ category_123update : chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
##  $ category_0321update: chr [1:869] "infectious" "infectious" "infectious" "infectious" ...
names(icd10categories)
## [1] "icd10code"           "category_original"   "category_1113update"
## [4] "category_123update"  "category_0321update"
dim(icd10categories)
## [1] 869   5
summary(icd10categories)
##   icd10code         category_original  category_1113update category_123update
##  Length:869         Length:869         Length:869          Length:869        
##  Class :character   Class :character   Class :character    Class :character  
##  Mode  :character   Mode  :character   Mode  :character    Mode  :character  
##  category_0321update
##  Length:869         
##  Class :character   
##  Mode  :character
sapply(icd10categories, class)
##           icd10code   category_original category_1113update  category_123update 
##         "character"         "character"         "character"         "character" 
## category_0321update 
##         "character"
glimpse(icd10categories)
## Rows: 869
## Columns: 5
## $ icd10code           <chr> "A00", "A01", "A02", "A03", "A04", "A05", "A06", "…
## $ category_original   <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_1113update <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_123update  <chr> "infectious", "infectious", "infectious", "infecti…
## $ category_0321update <chr> "infectious", "infectious", "infectious", "infecti…
colSums(is.na(icd10categories))
##           icd10code   category_original category_1113update  category_123update 
##                   0                 418                 418                   1 
## category_0321update 
##                   0
str(allcause_state_byweek)
## tibble [3,140 × 3] (S3: tbl_df/tbl/data.frame)
##  $ mmwr_week: chr [1:3140] "2019-01" "2019-01" "2019-01" "2019-01" ...
##  $ category : chr [1:3140] "cardiovascular" "digestive" "infectious" "infectious respiratory" ...
##  $ count    : chr [1:3140] "843" "441" "299" "240" ...
names(allcause_state_byweek)
## [1] "mmwr_week" "category"  "count"
dim(allcause_state_byweek)
## [1] 3140    3
summary(allcause_state_byweek)
##   mmwr_week           category            count          
##  Length:3140        Length:3140        Length:3140       
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character
sapply(allcause_state_byweek,class)
##   mmwr_week    category       count 
## "character" "character" "character"
glimpse(allcause_state_byweek)
## Rows: 3,140
## Columns: 3
## $ mmwr_week <chr> "2019-01", "2019-01", "2019-01", "2019-01", "2019-01", "2019…
## $ category  <chr> "cardiovascular", "digestive", "infectious", "infectious res…
## $ count     <chr> "843", "441", "299", "240", "691", "4090", "427", "220", "52…
colSums(is.na(allcause_state_byweek))
## mmwr_week  category     count 
##         0         0         0

Deficit Function Used

deficit_cumulative <- function(fit, start, end){
  if (!"curve_fit" %in% attr(fit, "type"))
    stop("This is not the correct deficit_model fit, needs curve fit.")

  ind             <- which(fit$date %in% seq(start, end, by = "day"))
  n               <- length(ind) # returns the # of elements in the ind object
  A               <- matrix(1, n, n)
  A[upper.tri(A)] <- 0                                         # It sets all the upper triangular elements of matrix A (excluding the diagonal) to 0
  A               <- sweep(A, 2, fit$expected[ind], FUN = "*") # This line multiplies each column (2 means columns) of the matrix A by a corresponding value from fit$expected[ind]
                                                               # So, now it's scaling each column of A by expected values from the model fit (fit$expected[ind])
  fhat <- matrix(fit$fitted[ind], ncol = 1)                    # creates a column vector (matrix with 1 column) called fhat from the values in fit$fitted[ind]

 
   # Below, this code is constructing a data frame with:
  
  fit_deficit   <- A %*% fhat                                    # multiply design matrix A by fitted values fhat to get cumulative modeled effect.
  obs_deficit   <- cumsum(fit$observed[ind] - fit$expected[ind])
  varcov       <- fit$x[ind,] %*% fit$betacov %*% t(fit$x[ind,]) # Compute variance-covariance matrix of the linear predictor based on model design matrix x and estimated betacov
  diag(varcov) <- fit$se[ind]^2                                  # Replace diagonal entries of varcov with squared standard errors — probably to correct for over/under-dispersion or model misspecification.
  fitted_se    <- sqrt(diag(A %*%  varcov %*% t(A)))
  sd           <- sqrt(diag(A %*% (fit$cov[ind,ind] + diag(fit$log_expected_se[ind]^2)) %*% t(A))) # this computes the variance of the observed cumulative deficit.
                                                                                                   # It adds fit$cov (covariance of observed data) and extra uncertainty from fit$log_expected_se.
  data.frame(date     = fit$date[ind], # The time points
             observed = obs_deficit,   # Cumulative observed - expected difference
             sd       = sd,            # Modeled cumulative deficit/excess from the model.
             fitted   = fit_deficit,   # Standard deviation of the observed cumulative deficit.
             se       = fitted_se)     # Standard error of the fitted cumulative deficit.
}

Deficit Model fits for the different categories

Cardio

exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 2, 20), 
                     by = "day")

# cardiovascular data
cardio <- all_cause_state_by_week %>%
  filter(category == "cardiovascular") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# deficit model with COVID-19 dates excluded
cardio_deficit <- compute_expected(cardio, 
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 8.85.
cardio_deficit
## # A tibble: 312 × 12
##    MMWRyear MMWRweek category       weekly_count_statelevel date      
##       <int>    <int> <chr>          <chr>                   <date>    
##  1     2019        2 cardiovascular 1233                    2019-01-06
##  2     2019        3 cardiovascular 1160                    2019-01-13
##  3     2019        4 cardiovascular 1201                    2019-01-20
##  4     2019        5 cardiovascular 1157                    2019-01-27
##  5     2019        6 cardiovascular 1286                    2019-02-03
##  6     2019        7 cardiovascular 1261                    2019-02-10
##  7     2019        8 cardiovascular 1141                    2019-02-17
##  8     2019        9 cardiovascular 1156                    2019-02-24
##  9     2019       10 cardiovascular 1208                    2019-03-03
## 10     2019       11 cardiovascular 1232                    2019-03-10
## # ℹ 302 more rows
## # ℹ 7 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>,
## #   population <dbl>, outcome <dbl>, log_expected_se <dbl>, expected <dbl>,
## #   excluded <lgl>
cardio_deficit <- cardio_deficit %>%
  mutate(date = as.Date(paste(MMWRyear, MMWRweek, 1, sep = "-"), "%Y-%U-%u"))
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `date = as.Date(paste(MMWRyear, MMWRweek, 1, sep = "-"),
##   "%Y-%U-%u")`.
## Caused by warning in `strptime()`:
## ! (0-based) yday 369 in year 2020 is invalid
combined_plot_1 <- ggplot(cardio_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Cardiovascular Hospitalizations in MA",
       x = "Date",
       y = "Outcome",
       color = "Data Point Status") +
  theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
combined_plot_1
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

fit_x <- cardio_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = as.Date("2019-01-14"),
               end = as.Date("2024-12-31"),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_x <- excess_plot(fit_x, 
            title = "Deviations from Expected Counts for Cardio in MA"
)
deficit_plot_x

fit_x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-09-23 2019-12-02       8.144139       8.669496    0.07650853    12062
## 2 2020-01-20 2020-02-03       8.276251       8.703193    0.14678711     3343
## 3 2020-03-02 2021-02-15       7.265876       8.874462    0.03594969    49893
## 4 2021-12-27 2022-01-31       8.304721       9.024090    0.10569035     6709
## 5 2023-05-01 2023-05-01       9.677494       9.573722    0.26665497     1303
##    expected     deficit        sd
## 1 12840.088   -778.0880 113.31411
## 2  3515.453   -172.4533  59.29126
## 3 60938.774 -11045.7736 246.85780
## 4  7290.145   -581.1446  85.38234
## 5  1289.028     13.9721  35.90303
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_1 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Cardio",
       title    = "Cumulative Deficit Hospitilization for Cardio Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_1

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_1_adj <- combined_plot_1 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_x_adj <- deficit_plot_x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_1_adj <- cumulative_deficit_plot_1 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_1 <- combined_plot_1_adj /
                       deficit_plot_x_adj /
                       cumulative_deficit_plot_1_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_1
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("Cardio_combined_figure.png",
       plot = final_combined_plot_1,
       width = 11,      
       height = 10.5,    
       dpi = 300,
       bg = "white")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("cardio_combined_figure.png",
       plot = final_combined_plot_1,
       width = 11,       # good for standard document width
       height = 10.5,    # leaves space for three stacked plots
       dpi = 300,
       bg = "white")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

Digestive

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 2, 27), 
                     by = "day")

# digestive data
digestive <- all_cause_state_by_week %>%
  filter(category == "digestive") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
digestive_deficit <- compute_expected(digestive, 
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 4.33.
digestive_deficit
## # A tibble: 312 × 12
##    MMWRyear MMWRweek category  weekly_count_statelevel date      
##       <int>    <int> <chr>     <chr>                   <date>    
##  1     2019        2 digestive 535                     2019-01-06
##  2     2019        3 digestive 570                     2019-01-13
##  3     2019        4 digestive 553                     2019-01-20
##  4     2019        5 digestive 576                     2019-01-27
##  5     2019        6 digestive 609                     2019-02-03
##  6     2019        7 digestive 599                     2019-02-10
##  7     2019        8 digestive 602                     2019-02-17
##  8     2019        9 digestive 562                     2019-02-24
##  9     2019       10 digestive 550                     2019-03-03
## 10     2019       11 digestive 622                     2019-03-10
## # ℹ 302 more rows
## # ℹ 7 more variables: statecount_numeric <dbl>, statecount_suppressed <lgl>,
## #   population <dbl>, outcome <dbl>, log_expected_se <dbl>, expected <dbl>,
## #   excluded <lgl>
# Plotting
combined_plot_2 <- ggplot(digestive_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly digestive ED Hospitilizations in MA",
       x = "Date",
       y = "Outcome",
       color = "Data Point Status") +
  theme_minimal()

combined_plot_2

fit_2x <- digestive_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = min(.$date),
               end = max(.$date),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_2x <- excess_plot(fit_2x, 
            title = "Deviations from Expected Counts for Digestive in MA"
)
deficit_plot_2x

fit_2x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2019-09-15 2019-10-13       3.986860       4.261767    0.07956446     2684
## 2 2020-03-01 2020-09-06       3.519113       4.442692    0.03432838    13267
## 3 2020-10-11 2021-03-14       3.714835       4.193817    0.03680019    11504
## 4 2021-12-05 2022-01-30       3.759757       4.186787    0.05877984     4556
## 5 2022-05-01 2022-05-22       4.368984       4.646495    0.09288423     2353
## 6 2022-09-18 2022-10-02       4.263148       4.543836    0.10606203     1722
##    expected    deficit        sd
## 1  2869.071  -185.0706  53.56371
## 2 16748.877 -3481.8769 129.41745
## 3 12987.299 -1483.2986 113.96183
## 4  5073.467  -517.4671  71.22827
## 5  2502.459  -149.4588  50.02458
## 6  1835.377  -113.3774  42.84131
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_2x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_2 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Digestive",
       title    = "Cumulative Deficit Hospitilization for Digestive Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_2

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_2_adj <- combined_plot_2 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_2x_adj <- deficit_plot_2x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_2_adj <- cumulative_deficit_plot_2 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_2 <- combined_plot_2_adj /
                       deficit_plot_2x_adj /
                       cumulative_deficit_plot_2_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_2

ggsave("Digestive_combined_figure.png",
       plot = final_combined_plot_2,
       width = 11,      
       height = 10.5,    
       dpi = 300,
       bg = "white")

Infectious

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 6, 5), 
                     by = "day")

# infectious data
infectious <- all_cause_state_by_week %>%
  filter(category == "infectious") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# deficit model with COVID-19 dates excluded
infectious_deficit <- compute_expected(infectious,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 2.66.
# Plotting
combined_plot_3 <- ggplot(infectious_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Infectious ED Hospitilizations in MA",
       x = "Date",
       y = "Outcome",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_3

fit_3x <- infectious_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = as.Date("2020-03-01"),
               end = as.Date("2024-12-31"),
               knots.per.year = 12,
               verbose = FALSE)

deficit_plot_3x <- excess_plot(fit_3x, 
            title = "Deviations from Expected Counts for Infectious in MA"
)
deficit_plot_3x

fit_3x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-05-03 2020-08-02       2.384625       2.718529    0.03797626     4495
## 2 2020-10-04 2021-05-30       2.231734       2.726452    0.02405327    10517
## 3 2022-01-23 2022-02-27       2.530161       2.785967    0.05872481     2044
## 4 2022-12-25 2023-01-22       2.565316       2.820316    0.06472516     1727
## 5 2023-06-18 2023-07-02       2.542539       2.771474    0.08283312     1027
## 6 2024-10-20 2024-12-15       2.525210       2.768028    0.04779399     3060
##    expected     deficit        sd
## 1  5124.405  -629.40467  71.58495
## 2 12848.349 -2331.34905 113.35056
## 3  2250.654  -206.65422  47.44106
## 4  1898.669  -171.66889  43.57372
## 5  1119.473   -92.47266  33.45852
## 6  3354.243  -294.24302  57.91583
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_3x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_3 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Infectious",
       title    = "Cumulative Deficit Hospitilization for Infectious Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_3

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_3_adj <- combined_plot_3 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_3x_adj <- deficit_plot_3x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_3_adj <- cumulative_deficit_plot_3 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_3 <- combined_plot_3_adj /
                       deficit_plot_3x_adj /
                       cumulative_deficit_plot_3_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_3

ggsave("Infectious_combined_figure.png",
       plot = final_combined_plot_2,
       width = 11,       
       height = 10.5,    
       dpi = 300,
       bg = "white")

infectious_respiratory

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 6, 12), 
                     by = "day")

# infectious respiratory data
infectious_respiratory <- all_cause_state_by_week %>%
  filter(category == "infectious respiratory") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
infectious_respiratory_deficit <- compute_expected(infectious_respiratory,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 0.699.
# Plotting
combined_plot_4 <- ggplot(infectious_respiratory_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Infectious Respiratory ED Hospitilizations in MA",
       x = "Date",
       y = "Outcome",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_4

# -- Fitting excess model
fit_4x <- infectious_respiratory_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = as.Date("2020-03-01"),
               end = as.Date("2025-12-31"),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_4x <- excess_plot(fit_4x, 
            title = "Deviations from Expected Counts for Infectious Respiratory in MA"
)
deficit_plot_4x

# -- Intervals of inordinate mortality found by the excess model
fit_4x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-03-29 2021-05-30      0.1241042      0.8568669    0.01013141     1036
## 2 2021-10-31 2022-03-20      0.4714432      1.3271874    0.02166535     1333
## 3 2023-01-15 2023-02-26      0.6535837      1.2568241    0.03651720      616
## 4 2023-04-09 2023-04-23      0.4208683      0.7014588    0.04167251      170
## 5 2024-11-03 2024-11-10      0.6424430      1.0037101    0.06105179      173
##    expected     deficit       sd
## 1 7152.9723 -6116.97228 84.57525
## 2 3752.6065 -2419.60646 61.25852
## 3 1184.5517  -568.55175 34.41732
## 4  283.3380  -113.33805 16.83265
## 5  270.2837   -97.28366 16.44031
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_4x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_4 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-infectious respiratory",
       title    = "Cumulative Deficit Hospitilization for infectious respiratory Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_4

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_4_adj <- combined_plot_4 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_4x_adj <- deficit_plot_4x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_4_adj <- cumulative_deficit_plot_4 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_4 <- combined_plot_4_adj /
                       deficit_plot_4x_adj /
                       cumulative_deficit_plot_4_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_4

ggsave("Infectious Respiratory_combined_figure.png",
       plot = final_combined_plot_4,
       width = 11,       # good for standard document width
       height = 10.5,    # leaves space for three stacked plots
       dpi = 300,
       bg = "white")

injury

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 4, 10), 
                     by = "day")

# injury data
injury <- all_cause_state_by_week %>%
  filter(category == "injury") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
injury_deficit <- compute_expected(injury,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 7.86.
# Plotting
combined_plot_5 <- ggplot(injury_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Injury ED Hospitilizations in MA",
       x = "Date",
       y = "Outcome",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_5

# -- Fitting excess model
fit_5x <- injury_deficit %>%                                                     
  deficit_model(exclude = exclude_dates,
               start = as.Date("2020-03-01"),
               end = as.Date("2024-12-31"),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_5x <- excess_plot(fit_5x, 
            title = "Deviations from Expected Counts for Injuries in MA"
)
deficit_plot_5x

# -- Intervals of inordinate mortality found by the excess model
fit_5x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-03-08 2020-09-20       6.347086       7.766319    0.04459825    24783
## 2 2020-10-04 2021-02-14       6.613078       7.488809    0.05273521    17808
## 3 2021-03-07 2021-04-11       7.120101       7.580370    0.09686766     5752
## 4 2021-05-09 2021-05-16       7.731598       7.976285    0.17210544     2082
## 5 2021-07-25 2021-08-08       7.979167       8.372155    0.14396843     3223
## 6 2021-12-19 2022-01-16       7.140402       7.682093    0.10682282     4807
## 7 2022-07-31 2022-08-21       8.268205       8.668967    0.12687117     4453
##    expected     deficit        sd
## 1 30324.572 -5541.57214 174.13952
## 2 20166.209 -2358.20891 142.00778
## 3  6123.830  -371.83007  78.25490
## 4  2147.891   -65.89052  46.34534
## 5  3381.738  -158.73836  58.15272
## 6  5171.673  -364.67274  71.91434
## 7  4668.838  -215.83814  68.32890
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_5x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_5 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Injury",
       title    = "Cumulative Deficit Hospitilization for Injury Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_5

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_5_adj <- combined_plot_5 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_5x_adj <- deficit_plot_5x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_5_adj <- cumulative_deficit_plot_5 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_5 <- combined_plot_5_adj /
                       deficit_plot_5x_adj /
                       cumulative_deficit_plot_5_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_5

ggsave("Injury_combined_figure.png",
       plot = final_combined_plot_5,
       width = 11,       # good for standard document width
       height = 10.5,    # leaves space for three stacked plots
       dpi = 300,
       bg = "white")

pneumonia

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2022, 5, 7), 
                     by = "day")

# pneumonia data
pneumonia <- all_cause_state_by_week %>%
  filter(category == "pneumonia") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
pneumonia_deficit <- compute_expected(pneumonia,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 2.32.
# Plotting
combined_plot_6 <- ggplot(pneumonia_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Pneumonia ED Hospitilizations in MA",
       x = "Date",
       y = "Outcome",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_6

# -- Fitting excess model
fit_6x <- pneumonia_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = as.Date("2020-03-01"),
               end = as.Date("2024-12-31"),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_6x <- excess_plot(fit_6x, 
            title = "Deviations from Expected Counts for Pneumonia in MA"
)
deficit_plot_6x

# -- Intervals of inordinate mortality found by the excess model
fit_6x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-04-26 2021-07-04       1.489661       2.618500    0.01756974    12636
## 2 2021-09-05 2022-05-22       1.947460       2.813849    0.02345135     9964
## 3 2022-07-17 2022-08-28       1.694437       2.049678    0.04663405     1597
## 4 2022-10-30 2022-11-13       2.562345       2.880032    0.08443981     1035
## 5 2023-01-22 2023-02-19       2.407862       2.978610    0.06651675     1621
## 6 2023-04-23 2023-05-21       2.294970       2.662542    0.06288869     1545
## 7 2024-01-28 2024-02-11       2.619286       2.974941    0.08581986     1058
##    expected    deficit        sd
## 1 22211.329 -9575.3287 149.03466
## 2 14396.797 -4432.7975 119.98666
## 3  1931.814  -334.8138  43.95240
## 4  1163.322  -128.3221  34.10751
## 5  2005.234  -384.2341  44.77984
## 6  1792.454  -247.4540  42.33738
## 7  1201.659  -143.6586  34.66495
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_6x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_6 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-pneumonia",
       title    = "Cumulative Deficit Hospitilization for pneumonia Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_6

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_6_adj <- combined_plot_6 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_6x_adj <- deficit_plot_6x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_6_adj <- cumulative_deficit_plot_6 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_6 <- combined_plot_6_adj /
                       deficit_plot_6x_adj /
                       cumulative_deficit_plot_6_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_6

ggsave("Pneumonia_combined_figure.png",
       plot = final_combined_plot_6,
       width = 11,       
       height = 10.5,    
       dpi = 300,
       bg = "white")

Renal

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 2, 27), 
                     by = "day")

# renal data
renal <- all_cause_state_by_week %>%
  filter(category == "renal") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
renal_deficit <- compute_expected(renal,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 2.42.
# Plotting
combined_plot_7 <- ggplot(renal_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Renal ED Hospitilizations in MA",
       x = "Date",
       y = "Outcome",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_7

# -- Fitting excess model
fit_7x <- renal_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = as.Date("2020-03-01"),
               end = as.Date("2024-12-31"),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_7x <- excess_plot(fit_7x, 
            title = "Deviations from Expected Counts for Renal in MA"
)
deficit_plot_7x

# -- Intervals of inordinate mortality found by the excess model
fit_7x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-03-08 2020-07-12       1.778983       2.317654    0.03009932     4551
## 2 2020-09-06 2020-09-27       2.194704       2.341201    0.06593235     1182
## 3 2020-11-01 2020-12-20       1.925472       2.186413    0.04505369     2074
## 4 2021-01-03 2021-02-28       1.965702       2.247830    0.04306949     2382
## 5 2022-01-02 2022-01-02       2.124147       2.335657    0.13170849      286
## 6 2022-05-15 2022-05-29       2.364289       2.607537    0.08034593      955
##    expected     deficit       sd
## 1 5929.0293 -1378.02930 77.00019
## 2 1260.8985   -78.89848 35.50913
## 3 2355.0697  -281.06966 48.52906
## 4 2723.8769  -341.87691 52.19077
## 5  314.4782   -28.47817 17.73353
## 6 1053.2544   -98.25443 32.45388
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_7x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_7 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-renal",
       title    = "Cumulative Deficit Hospitilization for renal Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_7

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_7_adj <- combined_plot_7 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_7x_adj <- deficit_plot_7x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_7_adj <- cumulative_deficit_plot_7 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_7 <- combined_plot_7_adj /
                       deficit_plot_7x_adj /
                       cumulative_deficit_plot_7_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_7

ggsave("Renal_combined_figure.png",
       plot = final_combined_plot_7,
       width = 11,       # good for standard document width
       height = 10.5,    # leaves space for three stacked plots
       dpi = 300,
       bg = "white")

Respiratory

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 6, 12), 
                     by = "day")

# respiratory data
respiratory <- all_cause_state_by_week %>%
  filter(category == "respiratory") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
respiratory_deficit <- compute_expected(respiratory,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 3.89.
# Plotting
combined_plot_8 <- ggplot(respiratory_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Respiratory ED Hospitilizations in MA",
       x = "Date",
       y = "Outcome",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_8

# -- Fitting excess model
fit_8x <- respiratory_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = as.Date("2020-03-01"),
               end = as.Date("2024-12-31"),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_8x <- excess_plot(fit_8x, 
            title = "Deviations from Expected Counts for Respiratory in MA"
)
deficit_plot_8x

# -- Intervals of inordinate mortality found by the excess model
fit_8x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-03-08 2021-05-30       2.675694       4.168126    0.02182343    23417
## 2 2021-11-21 2022-03-27       3.632627       4.441859    0.04166917     9293
## 3 2023-04-23 2023-05-07       4.186401       4.696742    0.10783183     1691
##   expected     deficit        sd
## 1 36478.39 -13061.3943 190.99318
## 2 11363.18  -2070.1786 106.59821
## 3  1897.14   -206.1403  43.55617
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_8x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_8 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-respiratory",
       title    = "Cumulative Deficit Hospitilization for respiratory Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_8

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_8_adj <- combined_plot_8 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_8x_adj <- deficit_plot_8x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_8_adj <- cumulative_deficit_plot_8 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_8 <- combined_plot_8_adj /
                       deficit_plot_8x_adj /
                       cumulative_deficit_plot_8_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_8

ggsave("Respiratory_combined_figure.png",
       plot = final_combined_plot_8,
       width = 11,      
       height = 10.5,   
       dpi = 300,
       bg = "white")

other

# COVID-19 exclusion dates
exclude_dates <- seq(lubridate::make_date(2020, 3, 1), 
                     lubridate::make_date(2021, 6, 5), 
                     by = "day")

# other data
other <- all_cause_state_by_week %>%
  filter(category == "other") %>%
    filter(!(MMWRyear == 2019 & MMWRweek == 1 | MMWRyear == 2025 & MMWRweek == 1)) %>%
  mutate(population = 7001399, # given from https://www.census.gov/quickfacts/fact/table/MA/PST045223
         outcome = statecount_numeric)

# excess model with COVID-19 dates excluded
other_deficit <- compute_expected(other,
                               exclude = exclude_dates,
                               frequency = 52)
## Overall death rate is 41.3.
# Plotting
combined_plot_9 <- ggplot(other_deficit, aes(x = date)) +
  geom_line(aes(y = expected), color = "blue", size = 1.5) +  # Expected outcome line
  geom_point(aes(y = outcome, color = factor(excluded)), size = 2) +  # Observed outcomes
  scale_color_manual(values = c("FALSE" = "green", "TRUE" = "red"),
                     labels = c("Included in Model", "Excluded from Model")) +
  labs(title = "Weekly Other ED Hospitilizations in MA",
       x = "Date",
       y = "Outcome",
       color = "Data Point Status") +
  theme_minimal()
combined_plot_9

# -- Fitting excess model
fit_9x <-  other_deficit %>% 
  deficit_model(exclude = exclude_dates,
               start = as.Date("2020-03-01"),
               end = as.Date("2024-12-31"),
               knots.per.year = 12,
               verbose = FALSE)

# -- Visualizing deviations from expected hospitalization in Massachusetts
deficit_plot_9x <- excess_plot(fit_9x, 
            title = "Deviations from Expected Counts for Other in MA"
)
deficit_plot_9x

# -- Intervals of inordinate mortality found by the excess model
fit_9x$detected_intervals
##        start        end obs_death_rate exp_death_rate sd_death_rate observed
## 1 2020-03-08 2021-05-23       35.63609       41.62600    0.06950269   307080
## 2 2021-08-22 2021-10-10       40.63174       42.79761    0.19933061    43766
## 3 2021-11-21 2022-02-13       36.00309       40.26895    0.15167805    63018
## 4 2022-03-06 2022-04-10       41.29956       43.11281    0.23101320    33364
## 5 2022-04-24 2022-05-08       42.59682       44.33594    0.33130392    17206
##    expected     deficit       sd
## 1 358695.65 -51615.6461 598.9121
## 2  46098.94  -2332.9402 214.7066
## 3  70484.75  -7466.7499 265.4896
## 4  34828.84  -1464.8444 186.6249
## 5  17908.48   -702.4761 133.8226
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_9x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_9 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-other",
       title    = "Cumulative Deficit Hospitilization for other Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_9

library(patchwork)

x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(from = x_limits[1], to = x_limits[2], by = "1 year")


combined_plot_9_adj <- combined_plot_9 +
  theme(legend.position = "bottom",
        legend.justification = c(1, 0),
        legend.box.just = "right") +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

deficit_plot_9x_adj <- deficit_plot_9x +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_deficit_plot_9_adj <- cumulative_deficit_plot_9 +
  labs(x = NULL) +
  scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y")

final_combined_plot_9 <- combined_plot_9_adj /
                       deficit_plot_9x_adj /
                       cumulative_deficit_plot_9_adj +
  plot_layout(ncol = 1, heights = c(1, 1, 1))


final_combined_plot_9

ggsave("Other_combined_figure.png",
       plot = final_combined_plot_5,
       width = 11,       
       height = 10.5,    
       dpi = 300,
       bg = "white")

Different Approuch

library(ggplot2)
library(patchwork)
library(cowplot)
library(scales)

# Common scale settings
x_limits <- as.Date(c("2019-01-01", "2025-01-01"))
x_breaks <- seq(x_limits[1], x_limits[2], by = "1 year")
category_labels <- c("Cardiovascular", "Digestive", "Infectious", "Infectious Respiratory",
                     "Injury", "Pneumonia", "Renal", "Respiratory", "Other")

# ----- CLEANERS -----
clean_combined <- function(plots, labels) {
  mapply(function(p, label) {
    p +
      labs(title = label, x = NULL) +
      scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
      theme_minimal() +
      theme(legend.position = "none",
            plot.title = element_text(hjust = 0.5, size = 11),
            axis.text.x = element_text(size = 8))
  }, plots, labels, SIMPLIFY = FALSE)
}

clean_deficit <- function(plot_list, labels) {
  mapply(function(p, label) {
    p +
      labs(title = label, x = NULL) +
      scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
      theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5, size = 11),
            axis.text.x = element_text(size = 8))
  }, plot_list, labels, SIMPLIFY = FALSE)
}

clean_cumulative <- function(plots, labels) {
  mapply(function(p, label) {
    p +
      labs(title = label, x = NULL, y = NULL, subtitle = NULL) +
      scale_x_date(limits = x_limits, breaks = x_breaks, date_labels = "%Y") +
      scale_y_continuous(labels = comma) +
      theme_minimal() +
      theme(plot.title = element_text(hjust = 0.5, size = 11),
            axis.text = element_text(size = 8))
  }, plots, labels, SIMPLIFY = FALSE)
}

# ----- APPLY CLEANING -----
combined_clean <- clean_combined(list(
  combined_plot_1, combined_plot_2, combined_plot_3,
  combined_plot_4, combined_plot_5, combined_plot_6,
  combined_plot_7, combined_plot_8, combined_plot_9
), category_labels)

deficit_clean <- clean_deficit(list(
  deficit_plot_x, deficit_plot_2x, deficit_plot_3x,
  deficit_plot_4x, deficit_plot_5x, deficit_plot_6x,
  deficit_plot_7x, deficit_plot_8x, deficit_plot_9x
), category_labels)
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
cumulative_clean <- clean_cumulative(list(
  cumulative_deficit_plot_1, cumulative_deficit_plot_2, cumulative_deficit_plot_3,
  cumulative_deficit_plot_4, cumulative_deficit_plot_5, cumulative_deficit_plot_6,
  cumulative_deficit_plot_7, cumulative_deficit_plot_8, cumulative_deficit_plot_9
), category_labels)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
# ----- LEGEND FOR COMBINED -----
legend_shared <- get_legend(
  combined_plot_1 +
    theme(legend.position = "bottom",
          legend.title = element_text(size = 10),
          legend.text = element_text(size = 9))
)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
## Warning in get_plot_component(plot, "guide-box"): Multiple components found;
## returning the first one. To return all, use `return_all = TRUE`.
# ----- FINAL COMBINED PLOT (balanced, readable) -----
combined_final <- wrap_plots(combined_clean, ncol = 3) +
  plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

combined_final <- combined_final +
  plot_annotation(
    title = "Weekly Observed vs Expected Hospitalizations by Cause",
    theme = theme(plot.title = element_text(hjust = 0.5, size = 15))
  )

ggsave("combined_all_categories.png",
       plot = combined_final,
       width = 12, height = 10, dpi = 300, bg = "white")
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
deficit_final <- wrap_plots(deficit_clean, ncol = 3) +
  plot_annotation(
    title = "Deviations from Expected Weekly Admissions",
    theme = theme(plot.title = element_text(hjust = 0.5, size = 14))
  )

ggsave("deficit_all_categories.png",
       plot = deficit_final,
       width = 11, height = 9.5, dpi = 300, bg = "white")


cumulative_grid <- wrap_plots(cumulative_clean, ncol = 3)

y_label <- ggdraw() +
  draw_label("Cumulative Deficit in Hospitalizations", angle = 90, size = 12)

cumulative_final <- plot_grid(
  y_label, cumulative_grid,
  rel_widths = c(0.04, 0.96),
  nrow = 1
)

cumulative_final <- cumulative_final +
  plot_annotation(
    title = "Cumulative Deficit in Hospital Admissions by Cause (2020–2024)",
    theme = theme(plot.title = element_text(hjust = 0.5, size = 15))
  )

ggsave("cumulative_all_categories.png",
       plot = cumulative_final,
       width = 12, height = 10, dpi = 300, bg = "white")

combined_final
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).

deficit_final

cumulative_final

ggsave("deficit_plot_x.png", plot = deficit_plot_x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_2x.png", plot = deficit_plot_2x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_3x.png", plot = deficit_plot_3x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_4x.png", plot = deficit_plot_4x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_5x.png", plot = deficit_plot_5x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_6x.png", plot = deficit_plot_6x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_7x.png", plot = deficit_plot_7x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_8x.png", plot = deficit_plot_8x, width = 8, height = 6, dpi = 300)
ggsave("deficit_plot_9x.png", plot = deficit_plot_9x, width = 8, height = 6, dpi = 300)

Cumulative deficit interactive plots

cumulative_deficit <- deficit_cumulative(fit_x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for Cardio Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-Cardio", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly
# -- Computing excess Hospitilization counts in Massachusetts from March 1, 2020 to May 9, 2020
cumulative_deficit  <- deficit_cumulative(fit_2x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

# -- Visualizing cumulative excess deaths in MA
cumulative_deficit_plot_2 <- cumulative_deficit %>%
  ggplot(aes(date)) +
  geom_ribbon(aes(ymin = observed- 2*sd, ymax = observed + 2*sd), alpha = 0.5) +
  geom_line(aes(y = observed),
            color = "white",
            size  = 1) +
  geom_line(aes(y = observed)) +
  geom_point(aes(y = observed)) +
  scale_y_continuous(labels = scales::comma) +
  labs(x        = "Date",
       y        = "Cumulative Deficit Hospitilization-Digestive",
       title    = "Cumulative Deficit Hospitilization for Digestive Cases in MA",
       subtitle = "During the first wave of Covid-19")
cumulative_deficit_plot_2

cumulative_deficit <- deficit_cumulative(fit_2x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for Digestive Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-Digestive", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#3

cumulative_deficit <- deficit_cumulative(fit_3x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for Infectious Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-Infectious", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#4

cumulative_deficit <- deficit_cumulative(fit_4x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for infectious respiratory Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-infectious respiratory", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#5

cumulative_deficit <- deficit_cumulative(fit_5x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for injury Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-injury", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#6

cumulative_deficit <- deficit_cumulative(fit_6x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for pneumonia Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-pneumonia", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#7

cumulative_deficit <- deficit_cumulative(fit_7x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for renal Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-renal", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#8

cumulative_deficit <- deficit_cumulative(fit_8x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for respiratory Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-respiratory", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly

#9

cumulative_deficit <- deficit_cumulative(fit_9x, 
                                        start = make_date(2020, 03, 01),
                                        end   = make_date(2024, 12, 28))

direct_plotly <- plot_ly(data = cumulative_deficit) %>%
  add_ribbons(x = ~date, 
              ymin = ~observed - 2*sd, 
              ymax = ~observed + 2*sd,
              name = "95% CI",
              fillcolor = "rgba(180, 180, 180, 0.5)",
              line = list(color = "transparent"),
              hoverinfo = "none") %>%

  add_lines(x = ~date, 
            y = ~observed,
            name = "Observed",
            line = list(color = "black", width = 2),
            hoverinfo = "text",
            text = ~paste("Date:", date, 
                         "<br>Cumulative Deficit:", round(observed, 1),
                         "<br>Lower CI:", round(observed - 2*sd, 1),
                         "<br>Upper CI:", round(observed + 2*sd, 1))) %>%

  add_markers(x = ~date, 
              y = ~observed,
              name = "Data Points",
              marker = list(color = "black", size = 4),
              hoverinfo = "text",
              text = ~paste("Date:", date, 
                           "<br>Cumulative Deficit:", round(observed, 1),
                           "<br>Lower CI:", round(observed - 2*sd, 1),
                           "<br>Upper CI:", round(observed + 2*sd, 1))) %>%
  layout(title = list(text = "Cumulative Deficit Hospitalization for other Cases in MA<br><sup>During the first wave of Covid-19</sup>"),
         xaxis = list(title = "Date"),
         yaxis = list(title = "Cumulative Deficit Hospitalization-other", 
                     tickformat = ","),
         hoverlabel = list(bgcolor = "white", 
                          font = list(family = "Arial", size = 12)),
         hovermode = "closest")

direct_plotly
ggsave("cumulative_deficit_plot_1.png", plot = cumulative_deficit_plot_1, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_2.png", plot = cumulative_deficit_plot_2, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_3.png", plot = cumulative_deficit_plot_3, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_4.png", plot = cumulative_deficit_plot_4, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_5.png", plot = cumulative_deficit_plot_5, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_6.png", plot = cumulative_deficit_plot_6, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_7.png", plot = cumulative_deficit_plot_7, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_8.png", plot = cumulative_deficit_plot_8, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_9.png", plot = cumulative_deficit_plot_9, width = 8, height = 6, dpi = 300)
ggsave("cumulative_deficit_plot_1.png", plot = cumulative_deficit_plot_1, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_2.png", plot = cumulative_deficit_plot_2, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_3.png", plot = cumulative_deficit_plot_3, width = 5, height = 3, dpi = 300)

ggsave("cumulative_deficit_plot_4.png", plot = cumulative_deficit_plot_4, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_5.png", plot = cumulative_deficit_plot_5, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_6.png", plot = cumulative_deficit_plot_6, width = 5, height = 3, dpi = 300)

ggsave("cumulative_deficit_plot_7.png", plot = cumulative_deficit_plot_7, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_8.png", plot = cumulative_deficit_plot_8, width = 5, height = 3, dpi = 300)
ggsave("cumulative_deficit_plot_9.png", plot = cumulative_deficit_plot_9, width = 5, height = 3, dpi = 300)
# new 'category' column added
cardiovascular <- fit_x$detected_intervals %>%
  mutate(category = "cardiovascular")
digestive <- fit_2x$detected_intervals %>%
  mutate(category = "digestive")
infectious <- fit_3x$detected_intervals %>%
  mutate(category = "infectious")
infectious_respiratory <- fit_4x$detected_intervals %>%
  mutate(category = "infectious respiratory")
injury <- fit_5x$detected_intervals %>%
  mutate(category = "injury")
pneumonia <- fit_6x$detected_intervals %>%
  mutate(category = "pneumonia")
renal <- fit_7x$detected_intervals %>%
  mutate(category = "renal")
respiratory <- fit_8x$detected_intervals %>%
  mutate(category = "respiratory")
other <- fit_9x$detected_intervals %>%
  mutate(category = "other")

combined_intervals <- bind_rows(
  cardiovascular,
  digestive,
  infectious,
  infectious_respiratory,
  injury,
  pneumonia,
  renal,
  respiratory,
  other
)

combined_intervals <- combined_intervals %>%
  select(category, everything())

combined_intervals <- combined_intervals %>%
  mutate(across(c(where(is.numeric)), ~round(., 1)))

kable_table <- combined_intervals %>%
  kable(
    format = "html",
    caption = "Detected Intervals by Disease Category",
    align = c("l", rep("c", 9)),
    col.names = c(
      "Category", "Start Date", "End Date", "Observed Count Rate", 
      "Expected Count Rate", "SD Count Rate", "Observed Counts", 
      "Expected Counts", "Deficit", "SD"
    ),
    row.names = FALSE
  )

kable_table
Detected Intervals by Disease Category
Category Start Date End Date Observed Count Rate Expected Count Rate SD Count Rate Observed Counts Expected Counts Deficit SD
cardiovascular 2019-09-23 2019-12-02 8.1 8.7 0.1 12062 12840.1 -778.1 113.3
cardiovascular 2020-01-20 2020-02-03 8.3 8.7 0.1 3343 3515.5 -172.5 59.3
cardiovascular 2020-03-02 2021-02-15 7.3 8.9 0.0 49893 60938.8 -11045.8 246.9
cardiovascular 2021-12-27 2022-01-31 8.3 9.0 0.1 6709 7290.1 -581.1 85.4
cardiovascular 2023-05-01 2023-05-01 9.7 9.6 0.3 1303 1289.0 14.0 35.9
digestive 2019-09-15 2019-10-13 4.0 4.3 0.1 2684 2869.1 -185.1 53.6
digestive 2020-03-01 2020-09-06 3.5 4.4 0.0 13267 16748.9 -3481.9 129.4
digestive 2020-10-11 2021-03-14 3.7 4.2 0.0 11504 12987.3 -1483.3 114.0
digestive 2021-12-05 2022-01-30 3.8 4.2 0.1 4556 5073.5 -517.5 71.2
digestive 2022-05-01 2022-05-22 4.4 4.6 0.1 2353 2502.5 -149.5 50.0
digestive 2022-09-18 2022-10-02 4.3 4.5 0.1 1722 1835.4 -113.4 42.8
infectious 2020-05-03 2020-08-02 2.4 2.7 0.0 4495 5124.4 -629.4 71.6
infectious 2020-10-04 2021-05-30 2.2 2.7 0.0 10517 12848.3 -2331.3 113.4
infectious 2022-01-23 2022-02-27 2.5 2.8 0.1 2044 2250.7 -206.7 47.4
infectious 2022-12-25 2023-01-22 2.6 2.8 0.1 1727 1898.7 -171.7 43.6
infectious 2023-06-18 2023-07-02 2.5 2.8 0.1 1027 1119.5 -92.5 33.5
infectious 2024-10-20 2024-12-15 2.5 2.8 0.0 3060 3354.2 -294.2 57.9
infectious respiratory 2020-03-29 2021-05-30 0.1 0.9 0.0 1036 7153.0 -6117.0 84.6
infectious respiratory 2021-10-31 2022-03-20 0.5 1.3 0.0 1333 3752.6 -2419.6 61.3
infectious respiratory 2023-01-15 2023-02-26 0.7 1.3 0.0 616 1184.6 -568.6 34.4
infectious respiratory 2023-04-09 2023-04-23 0.4 0.7 0.0 170 283.3 -113.3 16.8
infectious respiratory 2024-11-03 2024-11-10 0.6 1.0 0.1 173 270.3 -97.3 16.4
injury 2020-03-08 2020-09-20 6.3 7.8 0.0 24783 30324.6 -5541.6 174.1
injury 2020-10-04 2021-02-14 6.6 7.5 0.1 17808 20166.2 -2358.2 142.0
injury 2021-03-07 2021-04-11 7.1 7.6 0.1 5752 6123.8 -371.8 78.3
injury 2021-05-09 2021-05-16 7.7 8.0 0.2 2082 2147.9 -65.9 46.3
injury 2021-07-25 2021-08-08 8.0 8.4 0.1 3223 3381.7 -158.7 58.2
injury 2021-12-19 2022-01-16 7.1 7.7 0.1 4807 5171.7 -364.7 71.9
injury 2022-07-31 2022-08-21 8.3 8.7 0.1 4453 4668.8 -215.8 68.3
pneumonia 2020-04-26 2021-07-04 1.5 2.6 0.0 12636 22211.3 -9575.3 149.0
pneumonia 2021-09-05 2022-05-22 1.9 2.8 0.0 9964 14396.8 -4432.8 120.0
pneumonia 2022-07-17 2022-08-28 1.7 2.0 0.0 1597 1931.8 -334.8 44.0
pneumonia 2022-10-30 2022-11-13 2.6 2.9 0.1 1035 1163.3 -128.3 34.1
pneumonia 2023-01-22 2023-02-19 2.4 3.0 0.1 1621 2005.2 -384.2 44.8
pneumonia 2023-04-23 2023-05-21 2.3 2.7 0.1 1545 1792.5 -247.5 42.3
pneumonia 2024-01-28 2024-02-11 2.6 3.0 0.1 1058 1201.7 -143.7 34.7
renal 2020-03-08 2020-07-12 1.8 2.3 0.0 4551 5929.0 -1378.0 77.0
renal 2020-09-06 2020-09-27 2.2 2.3 0.1 1182 1260.9 -78.9 35.5
renal 2020-11-01 2020-12-20 1.9 2.2 0.0 2074 2355.1 -281.1 48.5
renal 2021-01-03 2021-02-28 2.0 2.2 0.0 2382 2723.9 -341.9 52.2
renal 2022-01-02 2022-01-02 2.1 2.3 0.1 286 314.5 -28.5 17.7
renal 2022-05-15 2022-05-29 2.4 2.6 0.1 955 1053.3 -98.3 32.5
respiratory 2020-03-08 2021-05-30 2.7 4.2 0.0 23417 36478.4 -13061.4 191.0
respiratory 2021-11-21 2022-03-27 3.6 4.4 0.0 9293 11363.2 -2070.2 106.6
respiratory 2023-04-23 2023-05-07 4.2 4.7 0.1 1691 1897.1 -206.1 43.6
other 2020-03-08 2021-05-23 35.6 41.6 0.1 307080 358695.6 -51615.6 598.9
other 2021-08-22 2021-10-10 40.6 42.8 0.2 43766 46098.9 -2332.9 214.7
other 2021-11-21 2022-02-13 36.0 40.3 0.2 63018 70484.7 -7466.7 265.5
other 2022-03-06 2022-04-10 41.3 43.1 0.2 33364 34828.8 -1464.8 186.6
other 2022-04-24 2022-05-08 42.6 44.3 0.3 17206 17908.5 -702.5 133.8
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(webshot2)

# Create styled HTML table
html_table <- combined_intervals %>%
  kable("html",
        escape = FALSE,
        caption = "Detected Intervals by Disease Category",
        align = c("l", rep("c", ncol(combined_intervals) - 1)),
        col.names = c(
          "Category", "Start Date", "End Date", "Observed Count Rate", 
          "Expected Count Rate", "SD Count Rate", "Observed Counts", 
          "Expected Counts", "Deficit", "SD"
        )) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                full_width = FALSE,
                position = "center")

# Save as HTML and then convert to PNG
save_kable(html_table, "combined_intervals_table.html")
webshot("combined_intervals_table.html", "combined_intervals_table.png", vwidth = 1200, vheight = 800)